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ABSTRACT 

Cosmological shocks are a critical part of large-scale structure formation, and are responsible for 
heating the intracluster medium in galaxy clusters. In addition, they are also capable of accelerating 
non-thermal electrons and protons. In this work, we focus on the acceleration of electrons at shock 
fronts, which is thought to be responsible for radio relics - extended radio features in the vicinity of 
merging galaxy clusters. By combining high resolution AMR/N-body cosmological simulations with 
an accurate shock finding algorithm and a model for electron acceleration, we calculate the expected 
synchrotron emission resulting from cosmological structure formation. We produce synthetic radio 
maps of a large sample of galaxy clusters and present luminosity functions and scaling relationships. 
With upcoming long wavelength radio telescopes, we expect to see an abundance of radio emission 
associated with merger shocks in the intracluster medium. By producing observationally motivated 
statistics, we provide predictions that can be compared with observations to further improve our 
understanding of magnetic fields and electron shock acceleration. 

Subject headings: cosmology: theory — hydrodynamics — methods: numerical — cosmic rays — 
radiation mechanisms: nonthermal 



1. INTRODUCTION 

The assembly history of galaxy clusters are wrought 
with violent mergers, high Mach- number flows, and ex- 
treme plasma physical interactions. Much of this results 
from a wide range of cosmological structure formation 
shocks dMiniati et al.||2001bl |Ryu et al.|2003| |Pfrommer| 
et al. 20061 IKang et al.||2007UHoelt et al.||2UU8| ISkill man 
et al ]|2008| jVazza et al.||2009| |Paul et al.||2011j ). These 
shocks, however, do more than simply heat the inflowing 
plasma. They a lso accelerate electrons and ions to rel - 
ati vistic speeds ( Bell 1978; Bland ford fc Ostriker|1978|). 



See Drury (1983), Blandford & Kichler| ( |1987p , and |Jones| 
fc EiTi son (1991) for reviews. These relativistic particles 
then act as signatures of merger and shock activity. The 
relativistic protons have radiative loss (e.g. collisions, 
pion decay, inverse Compton scattering) times compa- 
rable to the Hubble time, and therefore remain in the 
intracluste r medium and contribute to the total pressur e 
of the gas ( |Miniati et al.|[2001b[ |Pfrommer et al.| [2006). 
On the other hand, relativistic electrons have relatively 
short lifetimes, on the order of a few hundred million 
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years, and spend the remaining part of their life emit- 
ting synchrotron radiation as they gyrate about magnetic 
field lines. 

Relativistic protons in the intracluster medium are, in 
principle, most easily observed through their collisional 
interactions with thermal ions, leading to pion decay s 
that end in gamma-ray emission ( jPfrommer et al.|2 006 ) . 
This emission is being studied with the Fermi satellite, 
and while preliminary results hint at low levels of rela- 
tivistic ions, long integrations of ind ividual clusters are 
still forthcoming ( |Aleksic et al.||2010| ). 
Relativistic electrons have been more extensively stud- 
synchrotron 



ied in several galaxy clusters through their 
radiation (jRo ttgering et al.| |1997[ |Orru 
t al 



et al. 2007 

|Bonafede et al.||2009Mvan Weeren et al.||2009b[ ). These 
electrons are most frequently associat ed with objects 
called radio relics ( Ensslin et al.|[l998| ), which have ex- 
tended radio emission in the cluster exterior, are associ- 
ated with shocks, and have moderately polarized radio 
emission with spectral indices of a « 1 — 2 for surface 
brightness S ~ v~ a ). They are also, by definition, not 
associated with active galactic nuclei (AGN). This spec- 
tral shape most likely indic ates that these electrons wer e 
recently shock-accelerated (Blandford & Eichler 1987). 
Radio halos, on the other hand, have low polarization 
and usually follow the X-ray morphology in the centers of 
clusters. They are thought to be associated with turbu- 
lent acceleration and/or older p opulations of previou sly 
shock- accelerated electrons (e.g. |Brunetti et aT||2001 ). 

The origin of the shock-accelerated electrons is be- 
lieved to be primarily due to diffusive shock accelera - 
tion (DSA), as described in |Blandford fc Eichlerl ( p87| . 
This is a first-order Fermi mechanism in which electrons 
are accelerated by reflecting off magnetic field pertur- 
bations created by plasm a effects in shock waves. Re- 
cent numerical studies by |Spitkovsky ( 2008 ) have shown 
success in reproducing this mechanism through ab-initio 
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simulations using particle-in-cell (PIC) methods, though 
the incoming plasma flow was at a much higher ve- 
locity than the shocks discussed here. Studies of non- 
relativistic flow are ongoing because of the difficulty as- 
sociated wit h the range in r elevant timescales in such 
simulations ( |Spitko vsky 20081). 

If these electrons are capable of reaching high enough 
energies, they will emit synchrotron radiation in the pres- 
ence of magnetic fields. It is widely believed that the 
cluster surroundings are magnetized at relatively low 
field strengths, on the order of microgauss in the ICM 
(iFusco-Femiano et al.|[2001| |Reph ^li 'fc Gruber] |2003| 
Kephaeli et al. '2006; Ryi iet aT]|2008[ ), although mea- 
surements ma de by Faraday rota t ion indicate ~ 10 x 



larger values (Govoni et al. 2006 Clarke et al. 2001). 



Because the shock- accelerated electrons are expected to 
have a power-law distribution in energy, the synchrotron 
emission will also have a power-distribution in frequency. 
There have been a number of recent studies of these 
radio relics in inter ferometric observations of nearby 



galaxy clusters ( IRottgerin g et al.||l997l |Orru et al.||2007 
Bonafede et al.||2009| |van Weeren et al.||2UU9b| ). Ho 
ever, the total number ol clusters with known r elics is 
still fairly small (~ 22) ( |van Weeren et aL||2QQ9b| ). This 
is in part due to the low surface brightness ol these ex- 
tended sources. 

Studies of these objects comes at a critical point in 
time with a number of upcoming improvements in radio 
astronomy capabilities. The VLA is currently being up- 
graded to the "Expanded VLA" (EVLA), which will be 
roughly a factor of 10 better in terms of surface bright- 
ness sensitivity due primarily to an increase in band- 
width of up to 1 GHz at 1.5 GHz (jNapier|2006[ ). Several 
other telescopes will be coming online in the near future, 
such as LOFAR, possibly the SKA, and even lunar far- 
side l ow frequency arrays ( jBurns fc LUNAR Consortium 
2009) allo wing an unprecedente d view ol the synchrotron 
U niverse ( |Rudnick et al. [2009 ). Given the low surface 
brightness and spectrum of the emission, low frequency 
observations (~ 1 GHz and below) are most effective in 
providing information on the radio relics and halos. 

Here we set out to model these radio relics using high 
resolution, adaptive mesh refinement (AMR) cosmolog- 
ical simulations. These simulations include both dark 
matter and adiabatic baryonic physics, and allow us to 
model the shock acceleration of electrons and produce 
observationally-relevant radio luminosity functions and 
scaling relationships between cluster parameters such as 
synchrotron power, mass, and X-ray luminosity. By us- 
ing large simulation volumes, we are able to show statis- 
tics for thousands of objects, from which individual mor- 
phological and evolutionary analyses can be carried out. 

In Section [2j we introduce our simulation and galaxy 
cluster sample set and describe our shock finding method 
and synchrotron emission models. In Section [3J we study 
the global statistics of radio relics through projections 
and phase diagrams. Then, in Section [4j we describe the 
properties of individual halos, and use their statistics to 
generate radio luminosity scaling relationships and lumi- 
nosity functions. We end with Section [5] where we discuss 
the implications for future surveys, the limitations of our 
model, and future directions. 

2. METHODS 
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Fig. 1. — Mass function of halos in relic64 and relic200 at z = 0. 
Halos are found using HOP with a minimum number of 30 dark 
matter particles. Conservative estimates of the low mass cutoff are 
1O 12 M and 2x ~ 1O 13 M for relic64 and relic200, respectively. 
This corresponds to 512 and 320 particles. The discrepancy at low 
mass for relic200 is because of the lack of resolution of low-mass 
halos and poor force resolution at small scales. Dashed lines show 
the lower limit for the resolved halos for relic64 and relic 200 in 
black and bl ue, respectively. Also shown are the fits from |Warren| 
|et al.| ( [2006|> . Because of the similarity between the two cosmolo- 
gies, these nts differ by less than thickness of the line in this mass 
range. 

2.1. Enzo 

All simulations were run using the Enzo cosmology 



code (|Bryan~ fc Norm an|1997a|b| |Norman fc Bryan|1999 
|0'Shea~et al.||2UU5a|b| ). While a lull description can be 
lound in the cited papers, we will review the key aspects 
that are of importance to this work. 

Enzo uses block-structured ad aptive mesh refinement 



(AMR; iBerger & Colella M989) as a base upon which 
it couples an Eulerian hydro dynamic solver for the gas 
with an N-Body particle mesh (PM) solver ( |Efstathiou| 



et a l.|1985 [Hockney & Eastwood|1988 ) for the dark mat 
Users hi 



;he 



ter. Users have the choice ol solving the hydrodynamics 
with several m ethods including the piecew ise parabolic 
method (PPM; |Woodward fc Colella 1984) e xtende d for 
cosmological applications by [Bryan et~al.| (| 199 5[) and 
the ZEU S finite-difference method ( |Stone fc: Norman 
|1992a|b| ). In this work we utilize both methods, and 
restrict our studies to adiabatic gas physics. 

The AMR method that Enzo uses breaks the simula- 
tion domain into rectangular solid volumes called grids. 
These grids contain many computational elements called 
cells that set the resolution scale. Each grid exists on a 
level of refinement determined by the spatial resolution 
of its cells that ranges from — Imax, where l m ax is de- 
fined by the user. Once cells within a given grid satisfy 
the refinement criteria (based on overdensity, minimum 
resolution of the Jeans length, local gradients of hydrody- 
namical quantities, shocks, or cooling time), a new grid 
is created at the next higher level. In our simulations, we 
refine on overdensity of the gas and dark matter fields. 

2.2. Simulations 

We will focus on two simulations that both use N-body 
dynamics for the dark matter and adiabatic baryonic 
physics. The first simulation, hereby denoted as relic64, 
has a comoving volume of (64 /i _1 Mpc) 3 with 256 3 root- 
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grid cells and up to 6 levels of additional refinement. 
The AMR is done by inserting a higher-resolution region 
wherever a cell satisfies the refinement criteria. Here we 
require a gas or dark matter over density (S = p/p where p 
is the average density of gas or dark matter, respectively) 
of 8 to refine. Because refinement effectively splits a cell 
into 8 cells, this ensures that cells on each level have sim- 
ilar amounts of mass. This allows for a peak spatial res- 
olution of 3.9 /i _1 kpc (comoving). The simulation uses 
the ZE US hydro dynamic solver , with initial conditions 
from an |Eisenstein fc Hu| ( |1999[ ) power spectrum with a 
spectral index n s = 0.97. The cosmological parameters 
used are Q M = 0.268, Q B 0.0441, Q C dm = 0.2239 
, Q A = 0.732, h = #o/(100 km s^Mpc" 1 ) = 0.704, 
and erg = 0.82. The dark matter mass resolution is 
1.96 x 10 9 H~ 1 Mq. The simulation was started at a red- 
shift of z = 99 and run until z — 0, using approximately 
300,000 cpu-hours on the Texas Advanced Computing 
Center (TACC) Ranger supercomputer. 

We use a second simulation with a larger volume 
(200 hr 1 Mpc) 3 with more modest resolution, rclic200^ 
to capture a higher mass range for our simulated galaxy 
clusters. As before, it uses 256 3 root-grid cells and 
up to 5 levels of AMR. It has a peak resolution of 
24.4 /i _1 kpc (comoving) and a dark matter mass res- 
olution of 6.23 x 10 10 /i _1 M . The simulation uses a 
slightly different cosmology of n s = 0.96, VLm = 0.279, 
Q B = 0.046, ttcDM = 0.2239 , n A = 0.721, h = 
iV( 100 km s _1 Mpc _1 ) = 0.701, and a 8 = 0.817, consis - 
tent with WMAP Year-5 results ([Komatsu et al.||2008[ ). 



This simulation uses the PPM hydro dynamic solver, was 
also run on the TACC Ranger, and took approximately 
100, 000 cpu-hours to complete. 

The mass functions of the two simulations are shown 
in Figure [I] To calculate the mass function we begin 
by fi nding all the halos using the halo-finding algorithm 
HOP [ETsenstein fc Hu t (1998), implemented in an 
analysis and visualization system written in Python, de- 
signed for use w ith the adaptive m esh refinement codes 
including Enzo ( |Turk et aL 2011). This method finds 
halos by "hopping" from one dark matter particle to its 
most dense neighbor until a particle is its own highest 
density neighbor. All particles that find the same densest 
particle are then grouped into a single halo. The relic64 
simulation contains lO n M — 1O 14 M halos, but in our 
analysis we only consider objects with masses above 
10 12 M , corresponding to ~ 510 dark matter particles. 
This run was primarily designed to have superb resolu- 
tion capable of capturing the morphology and structure 
of the relics and shocks. The relic200 simulation contains 
5 x 1O 12 M - 8 x 10 14 M o halos. For this simulation, we 
consider only halos above 2 x 1O 13 M , corresponding to 
~ 320 dark matter particles. This simulation is designed 
to study the statistics of medium-sized clusters. While 
neither of these simulations capture the most massive 
clusters in the Universe (e.g. Coma), they provide insight 
to radio relic origins, structure, and evolution. Studies 
of very large volume simulations are reserved for future 
work. 

2.3. Shock Finding 



http://yt.enzotools.org 



To identify the shocks that ultimately accelerate the 
electrons that emit synchrotron radiation, we need an 
accurate shock identification algorithm. For this we use 
the te mperature- jump method given in |Skillman et al. 
( |2008[ ). Here we present an overview of the method for 
completeness. We use the Rankine-Hugoniot tempera- 
ture jump conditions to derive the Mach number: 



T2 _ (5M 2 -l)(M 2 + 3) 
Ti ~ 16M 2 



(1) 



where T2 and T\ are the post-shock (downstream) and 
pre-shock (upstream) temperatures, respectively. M. is 
the Mach number using the upstream (pre-shock) gas. 

A cell is determined to have a shock if it meets the 
following requirements: 

V-iT<0, VT-VS>0, T 2 >T U p 2 >Pu (2) 

where v is the velocity field, T is the temperature, p is 
the densit y, and S = T I p^~ x is the entropy. In our anal- 
ysis, as in |Skillman et al.| ( |2008| ), we have set a minimum 
preshock temperature of T = 10 4 K since the low-density 
gas in our cosmological simulations is assumed to be ion- 
ized (a reasonable assumption at z < 6). Therefore, any 
time the pre-shock temperature is lower than 10 4 K, the 
Mach number is calculated from the ratio of the post- 
shock temperature to 10 4 K. 

Once a shock is found, we identify the cell with the 
most negative flow divergence, choosing from a ray 
aligned with the temperature gradient. Therefore, even 
if several cells in a row qualify as a shock, only the "cen- 
ter" of the shock is marked as a shock, and the tem- 
perature jump is taken from the full jump across all the 
shocked cells. This relieves problems when shocks are 
spread out over several cells, especially when not aligned 
with the coordinate axes. In addition, this method has 
been implemented to run "on-the-fTy" in Enzo so that 
post-processing of the data is not needed. This method 
then saves the Mach number and quantities such as the 
pre-shock density and temperature directly along with 
the other hydrodynamical quantities during simulation 
output. The unique feature of this shock- finder is its 
ability to accurately identify off-axis shocks within AMR 
simulations and quantify their Mach number even if the 
shock is identified as being spread out across several cells. 

2.4. Synchrotron Emission 

In order to estimate the synchrotron emission from the 
shock waves, we follow the method of Hoe ft fc Brii ggen 
Here we summarize the main features of the 
The first assumption is that the electrons are 



([2007). 

mode 



accelerated to a power-law distribution that is related to 
the Mach number from diffusive shock acceleration the- 
ory. These accelerated electrons form an extension to 
the thermal, Maxwellian distribution that has a power- 
law form and exponential cutoff related to balancing the 
acceleration and cooling times of the electrons. 

The accelerated electrons then emit in the radio 
through synchrotron radiation. Since we are not per- 
forming magnetohydrodynamic simulations, we assume 
that the magnetic field is governed by flux freezing such 
that the magnetic field strength is related to the den- 
sity by B = 0.1pG( 1Q -4c m -3 ) 2 ^ 3 ? where n is the number 
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Fig. 2. — Projection of several quantities through the entire relic64 simulation volume at 



0. Shown are mass- weighted density 



[g cm 3 ] (upper left), mass- weighted temperature [K] (upper right), radio emission- weighted Mach number (lower left), and radio flux 
density [erg s -1 Hz -1 cm -2 ] (lower right). 



density. This is a reasonable assumption even in m erg- 
ing clus ters , as was found by iRoettiger et al~ (1999). In 
Section |4.4| and Appendix [A] we explore the effect of a 
modified magnetic field model. 

The total radio power from a shock wave of area A, 
frequency v bs, magnetic field B, electron acceleration ef- 
ficiency £ e , electron power-law index s (n e oc E~ s ), post- 
shock electron density n e and temperature T 2 is ( |Hoeft| 
|fc Bruggen||2QQ7|) 



dP(v ohs ) 
dv 



6.4 x 10 M erg s' 1 Hz 

£e / ^obs 



Mpc 2 10" 4 cm- 3 

)" S/2 x (^) 3/2 



0.05 y lAGHz 

(g/MG) 1+(s/2) 
(Bcmb/hG) 2 + (B/liG) 2 



7keV 
*(M)(.3) 



Note that the radiation spectral index is related to the 
electron spectral index by a = (s — l)/2. Bqmb is de- 



fined as the magnetic field corresponding to the energy 
density of the CMB. It has a value of B = 3.47/iG(l+z) 2 , 
and accounts for the inverse Compton emission that is 
simultaneously cooling the electrons along with their 
synchrotron emission. The final term, is a di- 

mensionless quantity that contains dependencies on the 
shock Mach number such that at \I/(2.5) ~ 10 -3 and ap- 
proaches 1 for M > 10. It can be thought of as a shape 
factor that, together with £, defines the acceleration ef- 
ficiency as a function of Mach number. Note that for 
calculations with z > 0, we modify v b s — » v b s x (1 + z) 
since we observe the redshifted emission. 

In all of our analys i s, we use £ e = 0.005, which was 
found by 
known re 



Hoeft et al. (2008) to match radio emission in 
ics to similar mass clusters in their simulation. 
While the true value of this parameter is quite uncertain 
from current observational and theoretical constraints, 
the relationship between it and the total radio power is 
linear. Therefore, if we underestimate the electron accel- 
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eration efficiency by a factor of 10, it will lead to a de- 
rived radio power that is low by a factor of 10, and thus 
all relationships between radio power and other quan- 
tities (e.g. mass, x-ray luminosity) simply need to be 
rescaled. While we have chosen the dimensio nless shape 
factor SP(M) from |Hoeft fc Bruggen) ( |2007[ ), there are 
still uncertainties in the efficiency of acceleration as a 
function of Mach number. However, exploration of the 
effects of these uncertainties are beyond the scope of this 
work. We also ignore the effects of re-accelerated 7 ~ 200 
electrons from radio galaxies. 

3. GLOBAL PROPERTIES OF RADIO RELICS 

3.1. Full Box Projections 

To begin our study of radio relics, we first performed 
simple projections of the radio emission through the en- 
tire simulation volume. An example is shown in Figure [2] 
along with projections of density, temperature and Mach 
number. For quantities such as density, temperature and 
Mach number in an AMR simulation, we choose to weight 



each cell by a secondary quantity since a simple average 
along the line of sight for each cell would bias the most 
highly refined regions because of their increased number 
of cells. Therefore, we choose to weight the density and 
temperature fields by cell mass, and the Mach number 
by the radio emission. This has the effect of pulling out 
the values of density and temperature from the densest 
regions, and the shocks that contribute the most to the 
radio emission. For radio and X-ray fields, we project 
the emissivities [energy /time/ volume] without a weight, 
leading to final values with units of [energy /time/area]. 
Mathematically, a weighted projection (here along the 
z-axis) is defined by: 



Pz(x,y) 



J w(x, y, z)v(x, y, z)dz 
J w{x, y, z)dz 



(4) 



where w(x,y,z) is the weight quantity at that location 
and v(x,y,z) is the value of the projected quantity. To 
evaluate this integral in our AMR setting, the integral 
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2D Phase Diagram of the relic radio emission in temperature-Mach number space for the relic200 (left) and relic64 (right) 
The temperature is that of the cell in the center of the shock. The color indicates the amount of 1.4 GHz radio emission 
_1 / (Mpc/ h) 3 1 (dlog A4 dlogT)] for a given value of temperature and Mach number. 



traverses the box along cells that are at the highest re- 
finement for a given point in space, and ignores cells that 
are covered by more highly refined regions. This, like the 
bulk of our analysis, is done using yt : detailed above. 

In Figures |2pl we see that in general the radio emis- 
sion traces out the large scale structure seen in the den- 
sity projection. Additionally, the emission is also highly 
correlated with the temperature structure. However, the 
correlation with Mach number is more interesting. In the 
projection of Mach number, we see shocks with strengths 
up to M ~ 10 — 100 throughout the volume in filaments 
and cluster edges, whereas the peak radio emission only 
shows up in small, curved arcs within clusters. At the 
location of these arcs, the value of the Mach number pro- 
jection drops to values between A4 ~ 3 — 10. This shows 
that the strongest shocks which are most likely external, 
accretion, shocks are not responsible for the bright ra- 
dio emission, and that it is instead t he interior shocks 
flRyu et a l. 2003; Skil lman et al.| [2008), as was found by 
|Hoelt et al.| pU(Jo| ) with moderate strengths, that shine 
in the radio. Tins can be understood by the fact that it is 
the mass flux of gas through shocks that is most impor- 
tant since that determines the number of electrons that 
can be accelerated. Therefore, while the Mach number 
is much lower for the interior shocks, the shock veloc- 
ity stays roughly constant while the pre-shock density is 
much higher, yielding more accelerated electrons. In the 
projection of Mach number, this results in the appear- 
ance of "veins" lining the interior of the filaments, "arcs" 
in the periphery of the clusters, and "holes" in the cen- 
ters of the clusters. While they are decrements in the 
projection of Mach number, they are the bright areas in 
the radio emission. 

The lack of strong emission in the accretion shocks sug- 
gests that having a hot, dense plasma is more important 
than the Mach mincer of the shock. This can be under- 

% oc A n e ZuT*/ 2 , Bl+S/2 



mce 



stood by Equation |3 S 
and in most cluster situations B^ MB 



2a + 

8/3 



> B 2 and s = 



1 = 3, W e have that ^ ex n e B 5 / 2 oc n e (ne /3 ) 5/2 oc 

n e ' . This implies that since the density in the accretion 
shocks is ~ 10 2 — 10 3 times lower than that in a merger 
shock, the power emitted will be down by a factor of 
« 2 x 10 3 — 10 5 . Therefore the features we see observa- 
tionally are more likely to be related to merger shocks 
than accretion shocks. 

3.2. Phase Diagrams 

The second method we use to study the bulk prop- 
erties of the radio emitting plasma is phase diagrams. 
Such diagrams are the equivalent of a two-dimensional 
histogram. Here we use them to study the gas properties 
of the radio emitting regions. 

The structure of these diagrams is as follows. For 
a given simulation output, we construct x and y-axis 
bins that are equally spaced logarithmically in two fields. 
Within each of these 2D bins, we integrate the total 
amount of a given quantity such as radio emission. This 
integrated value is normalized by the comoving volume 
of the simulation in order to give a comparable value 
between different physical size simulations. We have 
found three particularly insightful quantities to exam- 
ine in a range of permutions: temperature, overdensity, 
and Mach number. 

We have constructed one such phase diagram, seen in 
Figure |4j in which the x-axis is the Mach number, the y- 
axis is temperature, and the bins are colored by the total 
radio emission in that bin. The total integrated emission 
is normalized by the volume of the simulation and the 
size of the bins. As such, one reads this figure as "At 
Mach number x and temperature y, there is z amount of 
radio emission per comoving Mpc/h per AlogA4AlogT '." 
The utility of these diagrams is demonstrated in Figure 
[4j where it is immediately clear that the bulk of the radio 
emission in both simulations originates from hot gas with 
T = 10 6 - 5 x 10 7 K, and Mach number M = 3-10. This 
reinforces our earlier hypothesis that the radio features 
are generated from interior shocks associated with merg- 



7 




10" 3 10" 2 10 _1 10 10 1 10 2 10 3 10 4 iaV 3 io 2 io _1 10° 10 1 10 2 10 3 10 4 10 5 
Over dens ity Over dens ity 

Fig. 5. — 2D Phase Diagram of the relic radio emission in temperature-overdensity space for the relic200 (left) and relic64 (right) simulation. 
The color indicates the amount of 1.4 GHz radio emission [erg s" - 1 Hz~ 1 /(Mpc//i) 3 /(dlogS dlogT)] for a given value of temperature and 
overdensity. 



ing sub clusters that have low Mach numbers but high 
mass and energy flux due to the high relative density, 
and therefore shock velocity, of cluster cores. Second, it 
points out that shocks with M. = 20 — 100 have little 
role at z = in producing appreciable radio emission. In 
fact, their integrated luminosity is a factor of 500 — 1000 
less than their low-Mach number counterparts. 

At first glance, one also picks out a diagonal structure 
in the relic200 phase diagram that seems to be an upper 
limit on the temperature for a given Mach number. This 
is a very interesting feature that has a simple explana- 
tion. We calculate our Mach number using a minimum 
pre-shock temperature of 10 4 K. Now, while the gas at 
the location of the shock is not necessarily the pre- or 
post-shock temperature, it is bounded by those two val- 
ues. This is because the shock location is based on the 
cell with the most convergent flow, not the location of 
the pre- or post-shock gas. Because of this, if gas with 
pre-shock temperature T\ < 10 4 K is being accreted, 
the gas at the location of the shock will have a maxi- 
mum temperature of T max < ^'l^^ 10 4 K. This 
maximum coincides perfectly with the diagonal feature. 
Therefore, any gas below this line is likely pristine gas 
(that is, gas that has not been previously shocked) being 
accreted onto filaments or clusters for the first time. This 
gives us a proxy for the relative amount of accretion in 
a simulation. 

We can then use this diagnostic to study the role of 
accretion in the relic64 and relic200 simulations. While 
relic64 does have a small amount of accretion, it is far be- 
low that of relic200. This is because of the different mass 
clusters present in each of the two simulations; relic200 
has clusters that are up to an order of magnitude more 
massive than in relic64- Recalling that the accretion ra - 
dius r s = scales with mass (Bondi & Hoyle 1944), 
the clusters in relic200 are able to pull in and accrete 
more gas than those in relic64- This behavior will be 
studied as a function of redshift below. 



As we did in the temperature-Mach number phase 
space, we now examine the behavior in the temperature- 
overdensity plane in Figure |5j where overdensity is de- 
fined as p/p, where p is the mean matter density of the 
Universe, &Mpcrit- Here our earlier findings are rein- 
forced - the strongest emission is coming from the dens- 
est, hottest regions in the simulations. 

Both simulations exhibit the same general properties, 
though relic200 has hotter gas again due to the larger 
clusters. This suggests that not only will the most mas- 
sive clusters likely be associated with the strongest ra- 
dio emission, but that the strongest features will arise 
from merger shocks passing through the centers of clus- 
ter, which seems to be the case observationally. Below an 
overdensity of ~ 10 — 30, the radio emission greatly de- 
creases. This strongly disfavors the possibil ity of seeing 
cluster accretion shocks in agreement with |Hoeft et al. 
(2008). This is compounded by the fact that these accre- 
tion features are more diffuse and therefore have reduced 
surface brightness compared to the more compact merger 
shocks, making them difficult to study observationally. If 
we then compare relic64 to relic200, we see that there 
is a relative absence of accreting gas in the relic64 sim- 
ulation, reinforcing our earlier findings that relic64 has 
less accretion due to the smaller mass halos compared to 
relic200. In the following section we study this accretion 
as an evolutionary tracer in more depth. 

3.3. Radio Emission as a Proxy for Cluster Accretion 

We have seen at z = that there is a relatively small 
amount of radio emission coming from cluster and fila- 
ment accretion shocks. We now investigate whether or 
not this holds for earlier times. Figure [6] is an analog to 
our prior phase diagrams, but we now show the evolu- 
tion of this phase diagram for relic200 back to z = 2.0. 
We have normalized the emission by comoving volume 
in order to avoid confusion with the expansion of space. 
Even with this taken out, we see that there is a strong 
evolutionary trend in the origin of the radio emission. 
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Fig. 6. — 2d Phase Diagram of the relic radio emission in temperature-Mach number space for the relic200 simulation for varying redshifts. 
From top left to bottom right, we have z = 2.0, 1.5, 1.0, 0.5, 0.25, 0.00 



At z = 2, we see that the emission from accretion 
shocks is comparable to, if not above, that of the inte- 
rior merger shocks. This roughly translates to an equal 
amount of thermal energy being processed by mergers 
and accretion shocks at z = 2 since the efficiency of ac- 
celeration does not differ dramatically between the two 
cases. The line that we previously identified with ac- 
cretion shocks is now quite strong, and there is even a 
dominant population within the accretion regime, cen- 
tered around 10 6 K and Mach numbers of 40 — 60. 

By z = 1, the emission from accretion shocks has 
dropped by a factor of 10 while the merger shocks have 
increased by a factor of 2 — 3. Finally, by z = 0, nearly 
all of the radio emission due to accretion shocks has dis- 
appeared while the interior shocks have increased by a 
factor of 10 from z = 2. This drop in radio emission 
from accretion shocks coincides with the Universe begin- 
ning to accelerate due to 'dark energy' at z = 0.75 — 1.0. 
This process manifests itself by both depleting the voids 
leading to less mass to accrete, and decoupling of clusters 
from the cosmic expansion. Therefore, instead of growing 
in mass, and growing the radius of influence, dark energy 
dominates and pulls all the remaining matter away from 
the cluster faster than it can grow. We note here that 
this assumes that there is no evolution in the relative 
magnetic field strengths or acceleration efficiencies be- 
tween merger and accretion shocks. It may be the case 
that the magnetic field strength at accretion shock loca- 
tions is lower or higher at early times compared to the 
interior of proto-clusters. However, since we do not have 



conclusive evidence to this evolution, we have chosen to 
adopt the simplest model that assumes no evolution. 

This decrease in radio emission from accretion shocks 
is similar to the depletion of cosmic ray proto n acceler- 
at ion a t z < 1 from accretion shocks, found in Skillman 
al. (12008). This behavior is a novel perspective to 
view the effects of Dark Energy. If the radio relics from 
these accretion shocks are observable in the future, one 
should see a decrease in their frequency and power as z 
decreases. 



4. INDIVIDUAL OBJECT PROPERTIES 
4.1. Cluster Projections 

In this section, we take the opposite approach from 
the previous section and examine projections of individ- 
ual objects. We begin with our list of halos and make 
radial profiles that start at the density peak and con- 
tinue to the previously found r2oo, defined here as the 
radius where S = p/p = 200. We call the mass enclosed 
within this radius the virial mass, and also record sev- 
eral other quantities, such as X-ray luminosity and radio 
power within the virial radius. 

We begin by demonstrating the power of having a large 
sample of clusters in a single simulation by projecting 
only the 51 most massive clusters at z = along the x- 
axis in relic64 in Figure [7] The width and depth of each 
individual projection here is 4 Mpc/h. 

The most important result gleaned from these images 
is the morphological properties of cluster structure. If 
we first examine the gas density (top), we see that while 
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Fig. 7. — Density (top), temperature (middle) and 1.4GHz radio emission (bottom) for the 51 most massive halos in the relic64 simulation 
at z = 0. Mass decreases from top left to bottom right (2.5 x 10 14 — 2.0 x 1O 13 M0). Each individual image is 4 Mpc/h across. All images 
are projections down the x-axis. 



there is some amount of substructure, the density is cen- 
trally concentrated. Since X-ray emission closely follows 
the density distribution of the gas, this implies that the 
X-ray emission will be brightest in the centers of clusters. 
However, the radio emission (bottom) is brightest on the 
edges of the clusters and has very little correlation with 
the density structure. Instead, it more closely follows 
the temperature structure (middle). This is because the 
temperature is more strongly affected by shocks than the 
density (recall p<il P\ < 4 from shock jump conditions for 
7 = 5/3). Note, however, that the emission is still con- 
fined within high density regions inside the virial radius. 
Immediately from these images, we expect radio emission 
to be anti-coincident with t he X-ray emission, as is seen 



■ >i mi 

in existing relic examples flGiacintucci et al.||2 008 ; 



van 



Weeren et aLl|2009a| [Bonalede et al.l|2UU9| |(Jlarke & En 
sslin||2QQ6p . This behavior implies that shocks are more 
likely to appear in radio imaging than in X-ray surface 
brightness maps. 

Also visible in the radio emission are common features 
such as arcs and rings. These features are due to merg- 
ing subclusters as their bow shocks propagate through 
the ICM. These shapes are similar to what is seen in ob- 
served radio relics. This similarity supports our claim 
that the morphology of these objects is relate d to the 
locati o n of s hocks, as was originally suggested inlEnsslin 



et al. (1998). In a few rare situations (here in 2-3 clus- 
ters), these arcs appear in the very center of the cluster. 
Because the surrounding medium is both hot and quite 
dense in these cases, the radio emission is very str ong 



This agrees with our previous results from Section 3.2 



where we found the bulk of the emission at late times to 
be in the hot, dense phase of the gas. 

4.2. Radio Power - Mass Relationship 



From the projections of individual halos in Figure [7| 
we can see that there is a general trend for the more 
massive halos to have higher radio emission (note masses 
decrease from the top left to bottom right). We now 
want to quantify this scaling relationship by studying 
the radio luminosity-mass relationship for the halos in 
our simulations. We begin with the earlier list of halos 
and use the virial quantities of each halo. For each halo, 
we use their total mass and 1.4 GHz radio power (inte- 
grated out to r2oo) to populate Figure [8j Second, for the 
distribution of halos, we now determine the linear- least 
squares fit to log (Pi. 4GHz) — Alog(M2oo) + B for all ha- 
los with M 20 o > 1O 13 M and 2 x 1O 13 M for the relicM 
and re/ic200 simulations, respectively. We choose to only 
fit halos above this minimum mass because at smaller 
scales additional physics such as cooling not included in 
our simulations would possibly strongly affect the emis- 
sion. Additionally, we do not capture small mass halos 
that are likely moving through these small clusters pos- 
sibly creating a large fraction of the total radio emission. 
Because our simulation data does not have a measurable 
uncertainty for a given radio power, we have to use an 
alternate method of determining the error estimates of 
our parameters. We first find the best fit parameters 
using a uniform weighting. By calculating the residuals 
for each point from this best-fit relation, we estimate the 
uniform error for each point as the standard deviation of 
this residual. We then fit the data again using this error 
to obtain the uncertainty estimates in each parameter. 
The values of these parameters are shown in Table [II 

Perhaps the most interesting result from Figure~[8] is 
the normalization as a function of redshift. At z = 0, a 
1O 14 M mass halo emits ~ 10 20 erg/s/Hz, while a simi- 
lar mass halo at z = 2 emits ~ 10 22 erg/s/Hz. This am- 
plifies our hypothesis that the merger state of the halo 
is very important. The 1O 14 M halo at z = 2 is one 
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Fig. 8. — 1.4 GHz Radio Luminosity-Mass Relationship of halos in relic64 (left) and relic200 (right) for redshifts 2 (top), 1 (middle), 
and (bottom). Shown in black are each individual halo. The lines represents a best power-law fit to all halos with M v i r > 1O 13 M0 and 
2 x 1O 13 M0 for the relicQA and relic200 using a least-squares fitting routine. 



of the most massive objects at that time and has likely 
recently formed, whereas the same mass halo at z = 
is a fairly common object that was likely formed some 
time ago. Additionally, the probability of a merger with 
a 1 : 1 mass ratio is very low for the largest halos at a 
given time. For both simulations at multiple redshifts, 
the radio power is correlated with the mass of the halo. 
This dependence is expected since the radio emission is a 
function of the temperature, density and magnetic field 
strength of the halo, all of which scale positively with the 
mass. 



The second result of these radio luminosity-mass rela- 
tionships is the large scatter around the best fit. We see 
that there can be scatter of up to 2 — 3 orders of magni- 
tude for a given mass cluster. This suggests that while 
the radio emission is correlated with mass, the merger 
state of the halo plays a major role in determining the 
radio power. As can be seen from the projections of 
these halos in Figure [7| the most radio luminous objects 
have very disturbed morphology and are undergoing ma- 
jor mergers. 

In Figure [9| we show the evolution of the most massive 
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Fig. 9. — 1.4 GHz Radio emission overlaid on 0.2-12 keV X-ray emission for 20 snapshots during merger activity of the largest cluster in 
the relic64 simulation from z=0.23 to z=0.05, a time span of 2 .08 billion years. Length units are comoving. 
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Fig. 10. Time evolution of 1.4 GHz radio and 0.2-12 keV X- 
ray luminosities of the merger shown in Figure [9] normalized to 
the luminosities at z = 0.229. The first frames are a result of a 
prior merger in its last stages. 

cluster in relic64, and track the total radio and X-ray 
luminosity as a function of redshift. The X-ray emission 
is calculate d using outputs of the Cloudy code flFerlandl 
et al. 1 19981) where we have adapted the method oi |Smith| 
et al.[ ( |2(J08 ) for radiative cooling to calculate frequency- 
dependent emission. This yields an X-ray emission for 
given temperature and density of the gas, and is shown 
in gray. The radio emission, shown in color, has been 
masked such that all values below 10 -21 erg /s Hz~ 1 cm~ 2 



are transparent, allowing for a view of the X-ray data 
and masking out radio features that are too faint to be 
observed. To help follow the evolution of the total radio 
and X-ray luminosities from the cluster, we plot their 
relative luminosities with respect to their values at z = 
0.23 in Figure [10} 

As we follow the evolution from z = 0.23, we see the 
evolution of a major merger where the two cores pass 
through each other at z = 0.22. The smaller halo is 
moving from the center towards the upper-right. As the 
shockwave builds through z = 0.17, the radio emission 
closely follows the X-ray brightness jump, as we would 
expect. By z = 0.16, the radio emission from the initial 
shock has decreased dramatically. While the initial shock 
has disappeared in the radio, a secondary shock has been 
created that moves in from the hot ICM into the wake the 
merger left behind from the middle towards the lower- 
right. By z = 0.15 (frame 9), this is the most luminous 
feature in the radio. At this time, the image is 7-8 times 
brighter in the radio than it was at z = 0.16 (frame 8), 
illustrating how strongly radio emission depends on the 
merger state of the cluster. As the halo evolves further, 
additional smaller objects fall into the ICM, but don't 
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get nearly as bright as major merger. By z = 0.1 (frame 
14), the integrated radio luminosity has dropped back to 
pre-merger levels. 

During this merger, the X-ray luminosity also in- 
creases, but the total X-ray emission only increases by 
50%, which again illustrates the difference in the X-ray 
and radio emission mechanisms. A detailed analysis of 
cluster evolution and merger state is reserved for a later 
study. 

4.3. Radio Power-X-ray Relationship 

While the X-ray and radio emission may not be coinci- 
dent in projection, we expect to see a correlation between 
the total X-ray and radio luminosity since they both sam- 
ple hot, dense gas. We start our analysis from the results 
of our radial profiles and examine the 0.2 — 12keV X-ray 
and 1.4 GHz radio emission within r2oo for relic64 and 



relic200 in Figure 11 



We again use a method of linear-least squares regres- 
sion described above to obtain a scaling relationship be- 
tween the X-ray and radio luminosity, the results of 
which are shown in Table [I] Again we see that while 
there is a clear trend with X-ray luminosity, large scatter 
in the individual clusters can dominate the relationship. 
This scatter likely comes from two sources. First, a clus- 
ter that is relatively relaxed will have significant X-ray 
emission, whereas the lack of shocks in such a scenario 
will necessarily lead to zero radio emission in this model. 
Second, the fractional increase in X-ray luminosity across 
a shock front is much less than the radio emission because 
the X-ray primarily depends on the density of the gas at 
cluster temperatures, whereas the radio scales with den- 
sity and temperature. 

As a function of redshift, the scaling relationship be- 
tween radio power and X-ray luminosity evolves much 
like the radio power-mass relationship. Objects with the 
same X-ray luminosity at early times are more likely to 
have much higher radio power than their low-redshift 
counterparts. Additionally, the strength of the relation- 
ship increases at low redshift significantly due to the 
stronger correlation with larger mass halos. 

Even though our constraints on the fit parameters seem 
quite small, one can argue the due to the large scatter 
there should be a broad range of values capable of pro- 
ducing "acceptable" fits. To give a basic understanding 
of how our fit parameters can vary, for our two simu- 
lations at z = 0, we have shown the results of fitting 
using two cuts on our underlying data. The solid line is 
the result of fitting the points using all halos with amass 
greater than 1O 13 M and 2 x 1O 13 M for the relic6A 
and relic200 simulations, respectively. The solid line 
with hash marks shows the result of only fitting points 
with Lx > I0 44 erg/s. As one can see, the slope varies 
quite dramatically. Therefore when comparing to ob- 
servational constraints, the selection function of the ob- 
served/simulated clusters is very important. 

We can compare our derived scaling relationships with 
observatio nal estimates from known radio relics. IFer-l 
|etti| ( |2002[ ) found X-ray luminosities and lAGHz radio 
power tor 9 Abell clusters. If we fit their data using 
our same least-squares regression technique, we obtain 
Piaghz L x 0±0 ' 5 , agreeing quite well within the un- 
certainties in our z = simulation data. However, the 
normalization for the real relics are much higher than 



our simulated relics. We explicitly plot these clusters 
on Figure [TT] A very important point from this is that 
the observed relics land on the high end of both the X- 
ray luminosity and radio power. This demonstrates the 
selection effects coming in to play, as we have only ob- 
served the brightest objects as of yet. This suggests that 
deeper observations of radio-quiet clusters should lead to 
the dis covery of low power ra dio relics. As another con- 
straint, Cassano et al. (2006) study this relationship for 
giant radio halos and find P\aghz ^ L x 74±0 ' 21 . While 
these giant radio halos are thought to be from turbulent 
re-acceleration of electrons, their origin is likely linked 
to the same driving forces (i.e. mergers) as the radio 
relics. We note here that our X-ray emission is likely un- 
derestimated due to our lack of radiative physics in our 
simulations. However, properly modeling galaxy forma- 
tion, metal pollution, cooling, and thermal feedback is 
beyond the scope of this study. 

4.4. Luminosity Function 

Now that we have explored how the radio emission 
varies as a function of mass and X-ray luminosity, we ask 
the question: "How many radio relics do we expect at a 
given luminosity in the Universe?" Here we attempt to 
answer this question by constructing an observationally- 
motivated radio luminosity function. This is done by 
calculating the cumulative number of objects brighter 
than a given luminosity. We have done so in Figure 12 
and normalized the count rates by h 3 (l + z) 3 Gpc~ 3 . 

As we expect from our radio luminosity-mass relation- 
ship, the overall shape of the luminosity function is sim- 
ilar to the cluster mass function presented in Section 



4.2 



At first glance this luminosity function is not very 
encouraging for observational studies because even our 
most luminous objects are difficult or impossible to cap- 
ture with current radio telescopes. However, this is pri- 
marily a result of the mass range in our current simula- 
tions. The known radio relics are associated with massive 
clusters with M > 1O 15 M , and our largest clusters in 
relic200 only begin to reach the 1O 15 M mark. Obtain- 
ing spatial resolution needed to capture the relics in a 
volume large enough to capture these very rare clusters 
is computationally difficult. We address this by combin- 
ing best estimates of the halo mass function with our 
radio luminosity-mass relationship. 

We begin by calculating the best fits for the radio 
luminosity-mass relationships for both of our sim ula- 
tions. Next, we take fits from |Warren et aL ( |2006[ ) for 
the mass function at redshift and 1. We then convert 
the mass in the mass function to the expected radio lu- 
minosity from our fits in Table [T] The results of this 
fitting are shown in the left panel of Figure [12] Because 
of the scatter in the radio-mass luminosity function, we 
are able to place rough lower and upper limits on the 
luminosity function. This scatter will be constrained by 
future simulations that cover a larger mass scale of galaxy 
clusters. 

As we have done for our other results in Appendix A, 
we varied the magnetic field model to examine its ef- 
fects on the luminosity function. The first parameter we 
changed is the normalization of the magnetic field, Bo- 
In Figure [12} we show B = {0.01, 0.03, 0.1, 0.3, 1.0}/iG. 



Since the emitted power is roughly proportional to B 5 / 2 
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Fig. 11. — PiAGHz — Lx relationship for halos in relic64 (left) and relic200 (right) for z = 2, 1,0 from top to bottom. Both radio and 
X-ray emissivity are integrated out to the virial radius for each halo. A best fit line is found for halos with M v i r > 1O 13 M0 and 2 x 1O 13 M0 
for the relic64: and relic200 using a least-squares fitting routine. For z = 0, we show fits to our data using a minimum mass (solid) as 
above and minimum X-ray luminosity of 10 44 erg/s (solid + crosses). Also shown are observational data (stars) from |Feretti| ( |2002| along 
with a best fit (dashed). 



(see Section 3.1), as we increase Bo the luminosity func- 
tion shifts quite dramatically to larger luminosities. At 
low values of Bq the increase is close to the expected 
increase of £> 5 / 2 , while at higher values B approaches 
BcMB, reducing the effect of the increased local field 
strength. The second variation was in the scaling of the 
magnetic field with respect to the electron density. The 
line labeled "B-Flat" corresponds to B = Bo, whereas 
"B-Scale" denotes B oc Bon e . In both cases we set 
Bq = 0.1/iG. With a uniform magnetic field, we see that 



the number of high-luminosity objects decreases dramat- 
ically, while the number of low luminosity objects in- 
creases slightly. This is understandable given that the 
highest luminosity objects come from the most massive 
clusters, which have the highest densities. In this case, 
the density doesn't correspond to higher magnetic fields, 
and the radio luminosity is diminished with respect to 
the adiabatic scaling. Similarly, in the "B-Scale" case, 
the magnetic field strength is even higher in the dense 
parts of the largest clusters, leading to a shallower slope 
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TABLE 1 

Best fit parameters for radio power scaling relationship with halo mass and 0.2-12 keV X-ray emission. Fitting functions 
are 1O jB (M 2 oo/1O 13 M ) a and 10 B (L x / 10 43 erg / s) A for Mass and 0.5-12 keV X-ray luminosity, respectively. 
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i n the luminosity function. By comparison, Hoeft et al. 
(2008) also found a slope of —2/3 using the same model 
as the B = O.lfiG line in 12, adding verification to both 
results. 

To determine the number of clusters for a given survey 
area and redshift depth, we integrate the cosmological 
volume out to z = 0.5 for a given survey area dQ, 



Fig. 12. — (Left) 1.4 GHz Radio Luminosity Function for c lusters in relic64 an d relic200 for z = 0, 1, in units of inverse comoving 
(Gpc/h) 3 . Shown in dashed lines are extrapolations using the Warren et al.| ( |2006| mass function and our Pi.4gh z — Mass scaling from 
Table IT] (Right) The Radio Luminosity Function for clusters in relic t>4 as a function of magnetic field model. 

a function of the emitter's redshift. Therefore when de- 
riving the radio luminosity of an object at redshift z, we 
use v = 1.4 GHz (1 + z). Because of this, the emitted 
power is actually decreased since P\agh% v~ s I 2 where 
s ~ 2 for strong shocks. The power emitted in the clus- 
ter's frame is therefore substantially larger than what is 
The similar luminosity function is 



shown in Figure 12 



dV 
dzdVt 



(z) = 



c (l + z) 2 D 2 
H E(z) 



(5) 



with 



E 2 (z) = n mfi (i + z) 3 + n A , 



(6) 

where Da is the angular diameter distance. The result 
of this is that an all-sky survey out to z = 0.5 covers 
26.1 (Gpc/h) 3 . In combination with our estimates from 
the relic200 simulation in Figure [l2j we expect to find 
180 (conservative) to 1000 (optimistic) clusters with a 
total radio luminosity of 10 25 W/ Hz within this cosmo- 
logical volume. 

We also see from Figure 12 that the luminosity function 
of halos increases from z = 1 to z = 0. However, if plot- 
ted using proper volumes, the factor of 8 brings the two 
luminosity functions much closer together. Therefore the 
proper number density of radio relics seems to be fairly 
constant through cosmic time. This is an unexpected 
result, and encouraging for moderate redshift studies of 
radio relics. We note here that the frequency at which 
telescopes receive this synchrotron emission changes as 



therefore a product of the increased merger and accretion 
activity at higher redshift compared to that at z = 0. 

5. DISCUSSION 

5.1. Comparison To Previous Work 

In order to compare our results to previous shock stud- 
ies, we have calculated the kinetic energy flux through 
shocks, as shown in figure [13] The left two panels show 
the kinetic energy flux as amnction of Mach number in 
the relic64 and relic200 simulations, respectively. The 
right panels instead show the radio emission as a function 
of Mach number. While the black lines denote all tem- 
peratures, we also show the breakdown in terms of the 
pre-shock temperature. The kinetic energy flux results 



here c an be directly compared to F i gure 6 of Ryu et aIT| 
(|2003l), Figure 10 of |Skmman et al.| pQQ8] ) ? Figure 11 of 



Vazza et ah] (|2009|) , and c an also be compa r ed a fter unit 



conversions to Figure 6 of Pfrommer et al. 



J2006P . Even 
vary in size and adopted 



though these simulations 
cosmological parameters, the similarities in the kinetic 
energy flux processed by shocks is quite strong. This 
suggests that the underlying shock characteristics are 
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quite well understood even if the particular radio emis- 
sion models vary. 

In the right panels of [l3j we can see that there is a 
larger difference between the two simulations presented 
here in terms of the radio emission. This is likely due 
to the varying mass scales present in the simulations. 
We have also subsampled the relic200 simulation into 
random (64Mpc/h) 3 domains and found that a major 
contribution is confined to one of these subdomains. 

One of the earliest studies of shocks in a cosmological 
context is found in |Miniati et al.| ( |2QQQ| ) , where the au- 
thors found similar shock structure and kinetic energy 
flux trends as is seen in this study, though in a unigrid 
context. They, too, found that intermediate Mach num- 
ber shocks are responsible for processin g the majority 
of kin etic energy. In a pioneering work, Min iati et al.| 
( 2001a]) studied the injection and evolution of cosmic 
ray electrons. Using a framework to follow the cosmic 
ray distribution, they presented a radio power - core 
temperature relationship that shows strong similarity to 
what we have found with respect to cluster mass and 
X-ray luminosity, including a larger amount of scatter 
from cluster-to-cluster. Even though the resolution was 
modest compared to studies here, many of the primary 
characteristics of the radio emission are similar. 

Muc h of our work presen ted here can be compared with 
that of Hoeft et al. (2008). We use the same radio emis- 
sion model, but instead apply it to AMR simulations as 
opposed to smoothed particle hydrodynamic simulations. 
In particular, we can compare our radio relic luminosity 
function in Figure [l2| to their Figure 9. After accounting 
for the different normalization, we find that we have more 
objects at ~ 10 25 W/Hz. However, this result is from a 
small number of objects in our simulations and therefore 
future simulations with a larger sample of galaxy clusters 
are needed. 



In Pfrommer et al. (2008), the authors studied the ac- 
celeration and emission properties of cosmic ray electrons 
and protons in a smoothed particle hydrodynamics set- 
ting which focused on a set of high resolution galaxy 
clusters. Many of the same characteristics of galaxy clus- 
ter radio relic emission that we found in this study are 
consistent with their results. The morphology of the ra- 
dio relic emission is very similar to our results, though 
since they follow the electron population through time 
the emission is more diffuse compared to our simulated 
cluste rs. In another paper in the same series, |Pfrommer| 
( |2008| ) study the scaling relationship between the radio 
synchrotron, gamma-ray, and inverse Compton emission 
from the same set of galaxy clusters. Their results when 
fitting the scaling relationship between radio synchrotron 
emission and cluster mass give a significantly shallower 
slope of 1 — 1.5. However, due to the small number of 
clusters in their study, it is difficult to tell if there is a 
meaningful difference between their results and the ones 
presented here. In future work it would be useful to run 
a series of high resolution AMR simulations using our 
methods to compare to their results. 

5.2. Implications For Future Surveys 

Our study has shown that nearly every cluster has ra- 
dio emission and displays signs of radio relics at some 
stage in their evolution. When and where this radio 
emission occurs, however, is very sensitive to the merger 



and evolutionary state of the cluster. Current studies of 
radio relics have been confined to pointed observations 
of nearby, massive clusters, often based on strong X-ray 
emission. While this observational strategy does conform 
to our general results found in mass and X-ray scaling re- 
lationships, we have determined that not all X-ray lumi- 
nous or massive clusters have significant relic emission. 
It is instead heavily biased towards disturbed, merging 
clusters. Because the surface brightness of these relics 
is low due to their extended nature, large surveys with 
near-future telescopes are unlikely to yield serendipitous 
discoveries of cluster radio emission. Instead, the focus 
should be on deep, multiwavelength, large field-of-view 
observations with sensitivity to extended diffuse radio 
emission of disturbed X-ray clusters. 

Additionally, studies must include regions away from 
the peak X-ray emission. As was seen in Section [4] radio 
emission from shocks is generally brightest at the edges 
of clusters, surrounding the X-ray emission. This pre- 
scribes a fairly difficult observational roadmap, but the 
potential benefits include studying fundamental plasma 
physics phenomena such as in-situ shock electron acceler- 
ation and magnetic field structure. By combining statis- 
tical and morphological studies of these objects, we can 
readily compare them with the high-resolution hydrody- 
namical simulations presented here. 

Our lumin osity function, derived from the |Warren| 
[eTH] §006) mass function, suggests what can be ex- 
pected liTTuture observational studies. With surface 
brightness sensitivity improvements of a factor of 10, we 
can expect to see an increase of a factor of 10-100 in the 
number of clusters with radio relics. Alternatively, an all- 
sky survey out to z = 0.5 should result in the discovery 
of ~ 200 clusters above a luminosity of 10 25 W/Hz. 

5.3. Limitations of the Models 

There are several limitations to this study. First, 
we have not performed self-consistent MHD simulations, 
and instead adopted a scaling relationship between the 
post-shock density and magnetic field. This leads to the 
absence of any magnetic field configurations imprinting 
their structure on the radio relics. We plan to incorpo- 
rate MHD simulations in future work. However, since 
B ~ Be mb in the cluster environments we study here, 
the magnitude of our radio emission should not change 
dramatically. We are also not presenting a self consis- 
tent view of magnetic field generation, evolution, and 
non-linear interactions with shocks or particle accelera- 
tion, all of which are poorly understood in a cosmological 
context. 

Second, due to our relatively small box sizes, we have 
not captured the largest objects in the universe, which 
are likely to produce the brightest radio signatures. To 
account for this, we have provided an estimate using an 
analytic mass function combined with mass-luminosity 
scaling relationships. However, self-consistently captur- 
ing these very massive clusters is important, and will 
be explored in future work. Third, we also only fol- 
low recently accelerated electrons and ignore the aging of 
electron populations or re-accelerated electrons, which is 
likely to be important for radio halos and steep-spectrum 
objects. A more realistic model would follow the electron 
population "on-the-fly" and modify the acceleration effi- 
ciency as a function of pre-existing electron populations. 
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Fig. 13. — Left Panels: relicQA and relic200 kinetic energy flux processed by shocks at z = 0. Right Panels: 1.4 GHz radio emission. 



Given that we assume radio emission only comes from 
electrons that have just been accelerated, this means that 
the radio relic luminosities shown are lower limits. 

While our spatial and mass resolutions are quite good 
in the relicQA simulation, it is likely that incre asing the 
resolution w ould have an affect on our results. Skillman 
|et al.| ( |2QQ8[ ) found that in terms of the kinetic energy pro- 
cessed oy^hocks, a peak spatial resolution of 3.9kpc/h 
and a mass resolution of roughly 1O 9 M was approach- 
ing a converged result, though perfect convergence was 
not seen. In relic64 our spatial resolution matches this 
value, while we are a factor of two above this mass resolu- 
tion. For re/ic200, we are likely not capturing all of the 
kinetic energy flux in low Mach number shocks, which 
would suggest a higher spatial and mass resolution sim- 
ulation would lead to somewhat higher radio luminosity 
emission. This will likely always be the case since any in- 
crease in mass resolution will lead to a greater sampling 
of the mass function, allowing one to follow the merger 
assembly of galaxy clusters more accurately. This should 
increase the frequency of merger shocks, increasing overal 
radio emission. 

Finally, the electron acceleration efficiency in Equation 

ts poorly constrained at present. With additional radio 
servations, particularly using next generation low fre- 
quency radio telescopes, along with new PIC simulations, 
we may be able to calibrate £ e to more accurate values. 
This could be important in scaling the radio luminosity 
function in Figure [12] and estimating the number of radio 
relics expected in clusters and/or sky surveys. 

6. CONCLUSIONS AND FUTURE DIRECTIONS 

We have carried out high resolution AMR cosmological 
simulations using our accurate shock finding algorithm 
with a radio emission model for shock-accelerated elec- 
trons to examine the properties of radio relics in galaxy 
clusters. From this model, our main results are: 

• We have produced synthetic radio maps of the large 
scale structure and cluster environments, showing 
the variety of radio relic morphologies and loca- 
tions. 

• Through the use of 2D phase diagrams, we have 
found that while there is radio emission from both 
merger (internal) and accretion shocks, the emis- 
sion from the hot, dense intracluster medium asso- 
ciated with the merger shocks dominate the total 



emission. This balance is redshift-dependent, with 
accretion shocks being more important at high red- 
shift. 

• We have generated scaling relationships using over 
2000 simulated halos that give insight to how radio 
emission scales with mass and X-ray luminosity. 
These relationships evolve with redshift and there 
is a large scatter for individual halos as a result of 
merger state. 

• By studying the time evolution of a cluster under- 
going a merger, we find that the radio emission is 
highly dependent on the merger state, varying on 
time scales of a few hundred million years. 

• We have produced a synthetic radio luminosity 
function that gives observational predictions for the 
number of clusters with radio relics. This can be 
used to compare to future observed cluster luminos- 
ity functions and as a test of synchrotron emission 
models. 

In future studies we plan to examine the merger his- 
tory and morphology of these objects in greater detail. 
The redshift evolution of individual clusters is likely to 
be heavily correlated with merger state, and a statistical 
study of this relationship is vital to future observational 
studies. Finally, in order to correctly model these ra- 
dio relics, we need to self-consistent ly follow the electron 
population, taking into account effects of particle aging 
and re-acceleration. Future studies will also examine 
larger cosmological volumes and implement techniques 
such as light cones. 
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APPENDIX 
VARYING MAGNETIC FIELD MODELS 

Here we briefly examine the impact of an alternate magnetic field model. We will concentrate on two of our results, 
and leave further in-depth analysis to future work. In both cases, we present 6 alternate models. The first four 
concentrate on changing the reference magnetic field value, Bq = {0.01,0.03,0.3, 1.0}/iG. The second two models test 
the dependence of the magnetic field with respect to the electron number density. For these two cases, we chose 
simple parameterizations that simply give a feel for how this relationship affects our results. In one case, we keep the 
magnetic field constant throughout the domain at B = Bq = 0.1/iG. In the second case, we scale the magnetic field 
proportionally with the number density B = Bo(n e /n e; o), with Bo = 0.1/iG and n e) o = 10 -4 cm -3 . 

First we examined the fundamental impact on the phase of the gas responsible for the radio emission. In Figure [14] 
where we show the relative emission to our fiducial parameters, the primary effect of lowering Bq is to decrease the total 
emission as roughly £? 5 / 2 , as expected. If we instead increase Bq, because it starts to become comparable to Bqmb^ 
the emission enhancement is preferably increased in the high temperature, low-Mach number regime corresponding to 
merger shocks. In the case where we have a flat magnetic field, the emission in the high temperature, high density, 
low Mach number regions is greatly diminished. However, the accretion shock emission is increased by a factor of 
100-1000. If we let the magnetic field scale proportionally to number density, the opposite is true. Accretion shocks 
have weaker radio emission, whereas the merger shocks are more luminous. 

If we instead examine the effect of a changing magnetic field model on the radio luminosity-mass relationship, we 
again find a coherent picture, as shown in Figure 15 The key concept is that as the magnetic field approaches Bcmb, 
the added emission per unit increase is diminished by the second B 2 term in the denominator of Eq. Therefore, 
for low values of Bq when the magnetic field is much lower than Bqmb^ the higher magnetic field strengths in larger 
clusters has more of an effect, steepening the radio luminosity-mass relationship. On the other hand, when Bq is 
increased, the relative gains in magnetic field in large clusters does not impact the emission as strongly, flattening the 
relationship. In changing the B oc n e relationship, a flat magnetic field removes the density bias between large and 
small clusters, thereby flattening the scaling relationship. A linear scaling with number density increases this bias, 
steepening the relationship. 
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Fig. 14. — Relative radio emission with respect to our fiducial magnetic field model for the relic200 simulation at z = 0, shown in Figure 
In the upper four panels, the reference magnetic field parameter, Bo is varied. In the bottom left panel, the magnetic field strength is 
nat at B = 0.1/iG. In the bottom right panel, the magnetic field strength scales linearly with electron number density. 
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Fig. 15. — As in Figure [8] but for varying magnetic field model parameters. In the upper four panels, the reference magnetic field 
parameter, Bo is varied. In the bottom left panel, the magnetic field strength is flat at B = 0.1/iG. In the bottom right panel, the magnetic 
field strength scales linearly with electron number density. 
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